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The distribution of Yang-Lee zeros in the ferromagnetic Ising model in both two and three di¬ 
mensions is studied on the complex field plane directly in the thermodynamic limit via the tensor 
network methods. The partition function is represented as a contraction of a tensor network and is 
efficiently evaluated with an iterative tensor renormalization scheme. The free-energy density and 
the magnetization are computed on the complex field plane. Via the discontinuity of the magneti¬ 
zation, the density of the Yang-Lee zeros is obtained to lie on the unit circle, consistent with the 
Lee-Yang circle theorem. Distinct features are observed at different temperatures—below, above 
and at the critical temperature. Application of the tensor-network approach is also made to the 
g-state Potts models in both two and three dimensions and a previous debate on whether, in the 
thermodynamic limit, the Yang-Lee zeros lie on a unit circle for q > 2 is resolved: they clearly do not 
lie on a unit circle except at the zero temperature. For the Potts models {q = 3,4, 5, 6) investigated 
in two dimensions, as the temperature is lowered the radius of the zeros at a fixed angle from the 
real axis shrinks exponentially towards unity with the inverse temperature. 
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I. INTRODUCTION 

More than half a century ago Yang and Lee proposed 
a new approach for studying phase transitions in a gas 
by examining the zeros of the grand partition function 
in the complex fugacity plane^. In the thermodynamic 
limit, these zeros, which shall be referred to as the Yang- 
Lee zeros, may lie arbitrarily close to certain points on 
the real axis, marking where phase transitions appear. 
Inside a region clear of zeros, no phase transitions can 
occur. Thus, the study of the location of the zeros in 
the complex plane determines the transition points in 
the real axis^. In a subsequent paper^, Lee and Yang 
showed, among other things, that the zeros of the parti¬ 
tion for the ferromagnetic Ising model he on a unit circle 
of the complex field (or more precisely the complex 2 - 
plane, where z = exp(—2/3h) with ff = l/keT and h 
the external field), which is referred to as the Lee-Yang 
theorem. The equation of the state can be obtained via 
the knowledge of the distribution of the Yang-Lee zeros 
on the unit circle. At very high temperatures, the zeros 
do not cover the whole circle, but only the segment of 
the arc around the angle 8 = tt. As the temperature is 
lowered to the transition temperature Tc, the zeros move 
and pinch the real axis at 0 = 0, the field value (in this 
case zero) of which corresponds to the phase transition 
point. 

The study of partition function zeros, as well as gen¬ 
eralization of the circle theorem, has been extended to 
higher spins and other models (whose list is hard to ex¬ 
haust here)'^“^® and has provided useful insights to phase 
transitions, even in QCD^'*. The consideration of zeros in 
the complex temperature plane was initiated by Fisher^^ 
and these zeros are called Fisher zeros^’^^’^^. Proper¬ 
ties of the distributions of zeros also give characterization 
of first-order phase transitions^^ as well as higher-order 
phase transitions and scaling relations^®. The Lee-Yang 


theorem has its impact also beyond statistical physics; 
it has incited mathematical theory connected with the 
Laguerre-Polya-Schur theory of linear operators, possi¬ 
ble connection to the Riemann hypothesis on the zeta 
functions^*^’^ ' , and computational complexity of comput¬ 
ing averages^®, etc. Even though the Yang-Lee zeros lie 
on the complex plane, their density was inferred from 
the magnetization data in one experiment^®, and very 
recently it is found that the individual zeros of classical 
Ising models can be detected by using a quantum spin as 
a probe that couples to the classical spins®®’®b 

Yang-Lee zeros can be solved exactly on small system 
sizes and the thermodynamic limit is inferred from ex¬ 
trapolation. For most models, however, there is no an¬ 
alytic expression for the distribution of the zeros and to 
locate these zeros accurately in the thermodynamic limit 
remains a challenge. It was known that computing av¬ 
erages such as the magnetization exactly for even the 
ferromagnetic Ising model on the complex plane is ^P 
hard^®. Here we introduce the tensor network (TN) tech¬ 
niques for computing the density of the Yang-Lee zeros 
in the complex field plane directly in the thermodynamic 
limit. The density is proportional to the discontinuity 
of the magnetization on the complex plane®, and the 
magnetization is computed with controlled approxima¬ 
tions. Tensor network methods, including the density 
matrix renormalization group (DMRG), matrix product 
states (MPS), and other tensor product or tensor net¬ 
work states, have become very useful numerical tools in 
both classical and quantum spin systems®^^®' . In essence, 
physical observables are expressed as contraction of a ten¬ 
sor network. However, the contraction in two and higher 
dimensions is a computationally hard problem, but it can 
be approximated and the error can be controlled by re¬ 
ducing truncation, such as during the real-space coarse- 
graining or renormalization group (RG) procedure over 
the tensors describing the system®®""*® . We note that 
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the particular partition function zero closest to the pos¬ 
itive real axis has been explored with TN for the one¬ 
dimensional Schwinger model on a finite lattice"^^. But 
for our purpose we use the method appropriate for the 
infinite system to directly probe the density of zeros in 
the thermodynamic limit. 

In particular, we study ferromagnetic Ising and q-state 
Potts models^® in both two and three dimensions. We 
observe clearly that the magnetization has a dicontinu¬ 
ity at the unit circle, consistent with the Lee-Yang circle 
theorem for the Ising model. The distribution of the ze¬ 
ros on the unit circle shows distinct features at different 
temperature regimes: T > Tc, T = Tc and T < Tc- At 
T > Tc the zeros occupy part of the circle and move along 
it towards the real axis as T decreases. They pinch z = 1 
at T = Tc and the difference with T < Tc is manifest in 
the distribution of the zeros. From the density we ob¬ 
tain good estimates for the magnetization-field exponent 
S for both two and three dimensions. For the q-state 
Potts models in two and three dimensions our results 
demonstrate that the zeros clearly do not lie on a unit 
circle in the thermodynamic limit for q > 2 and T > 0, in 
agreement with Kim and Creswick^®, whose results were 
questioned previously due to the small system sizes'^®. 
Our results also show that the deviation from the unit 
circle, in terms of the radius at a fixed angle, is decaying 
exponentially with the inverse temperature /3, consistent 
with the fact that the location of the zeros shrinks to the 
unit circle at zero temperature. 

The remaining organization of the paper is as follows. 
In Sec. II, we give a slightly more expanded but elemen¬ 
tary introduction to the Yang-Lee zeros, and the readers 
can skip it if they are already familiar with these ideas. 
In Sec. Ill we review the RG algorithm that we employ 
in this paper, with a detailed exposition of the thermo¬ 
dynamic limit calculations of the free energy density. A 
numerical RG process is applied to the Ising model in 
both a square and cubic lattices, with computation of the 
magnetization described in Sec. IV, to show the Yang-Lee 
zeros on the complex plane. From the magnetization we 
directly obtain the density of zeros along the unit circle. 
In particular the distribution of the zeros at T = Tc is 
used to estimate the magnetization-field exponent <5. For 
T > Tc we study in Sec. V the ‘edge singularity’ of the 
density of zeros. In Sec. VI we extend our results to the 
Potts models, and we study the location of zeros in 2D 
and 3D lattices. We conclude in Sec. VII. 


II. YANG-LEE ZEROS: AN ELEMENTARY 
INTRODUCTION 


where Si = ±1 is a classical variable and h is an external 
field. The partition function at a temperature T = 1//3 
(where we have set fcs =1) is 


N 

Z(/3,/i)=Tr(e-^") = e'^^'‘^P„2", (2) 

n=0 


where in the second equality we have re-arranged the con¬ 
tribution to Z in terms of the number n of down spins 
and their associated value Pn at zero field and we have 
defined z = exp{—2j3h). Since Pn is independent of h 
(and is real and positive), one can factorize the polyno¬ 
mial V{z) = J2n=o PriZ"" = Co “ ^ri), where z„’s 

are the zeros of P{z), which are referred to as the Yang- 
Lee zeros and cq is a positive constant. 

Lee and Yang proved that for any ferromagnetic Ising 
model the zeros lie on a unit circle, namely = e*®”. 
This is referred to as the Lee-Yang circle theorem^. The 
consideration of the zeros and the Lee-Yang theorem was 
generalized to other models and higher spins'^”^^. 

How are the Yang-Lee zeros related to phase transi¬ 
tions? The free-energy density is 
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In a region free of zeros, the free energy is analytic and 
therefore there cannot be any singularity, and hence on 
the real axis contained in this zero-free region, no phase 
transitions occur. The Yang-Lee zeros close to the posi¬ 
tive real axis of z signal the location of the phase transi¬ 
tions. From the free energy, one can obtain the magneti¬ 
zation (which in general has a complex value for complex 

z), 


n—1 

In the limit ?V —> oo, the zeros form a continuum on the 
unit circle, with the density of zero per angle denoted 
by g{0). Since any complex zero Zn has a corresponding 
complex conjugate partner zjf, the density of zero has the 
symmetry that g{—d) = g{0)^ and, employing this, one 
can re-rewrite the free-energy density as follows, 

f{f3,h) = —h—^ f d9 g{9)\n{z^ — 2zcos9 + 1),{6) 
P Jo 

and the magnetization as follows 


Here we give an elementary introduction to the Yang- 
Lee zeros. Readers who are familiar with Yang-Lee zeros 
can skip this section. If we consider a system of N Ising 
spins with the Hamiltonian, 

if = - ^ SjSj - Si, (1) 

(id) ® 


= 1 — Az f d9 g{9)- 
Jo ' 


z — cos ( 


— 2zcos0 -I-1’ 
which we have chosen the normalization that 


n2-K 

/ d9g{9) = l. 
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( 7 ) 

( 8 ) 
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The thermodynamics of the system therefore can be de¬ 
duced if the density of zeros g{9) is known on the unit cir¬ 
cle. We remark that, however, except in one dimension, 
the analytic expression of g{9) is generally not known. 
But the Yang-Lee picture for phase transitions is very 
visual and intuitive in terms of partition function zeros; 
see e.g. Figs. 3, 5& 6. Approximation schemes such as 
series expansion have been applied to the density®^. Us¬ 
ing an electrostatic analogy, Lee and Yang showed that 
the density of zeros is related to the discontinuity of the 
magnetization across the unit circle, i.e., 

g(9) = -Ref lim miz = re^^)— lim m(z = re^^)] 

47r r^l— J 

(9) 

Since the free-energy density is expressed in terms of the 
zero density, other thermodynamic quanities in addition 
to the magnetization and the susceptibility, such as the 
specific heat and the entropy density can also be obtained 
once the density is known. Our numerical calculation, 
to be discussed below, is based on the above relation be¬ 
tween the magnetization and the density and is to obtain 
the latter via evaluation of the magnetization directly 
on the complex plane with the tensor network methods. 
We remark that we shall loosely refer to the real part of 
the magnetization as the magnetization when no confu¬ 
sion arises, as the imaginary part of the magnetization is 
never needed in our discussions. 

More than a decade after the results by Lee and Yang ’, 
Fisher initiated the study of the zeros defined on the com¬ 
plex temperature plane, i.e., the so-called Fisher zeros’^L 
Analyzing the behavior of Yang-Lee and Fisher zeros also 
allows characterization of first-order phase transitions^"* 
as well as second and higher-order phase transitions and 
scaling relations^®. The interest of Lee-Yang-like theo¬ 
rem for partition function zeros goes beyond statistical 
physics, such as in mathematics'^®’^^ and computational 
complexity'^®. But further detailed explanation on these 
will take us awry from the main purpose of this work 
and the readers are referred to the cited references and 
more references therein for further discussions. This sec¬ 
tion serves as a brief and elementary introduction to the 
Lee-Yang zeros that we shall explore in later sections. 


III. FREE ENERGY FROM TENSOR 
RENORMALIZATION 

The partition function Z = Tr exp{—/3iL} of a classi¬ 
cal system with local interactions can be expressed as a 
contraction of a tensor network of low rank and of low 
dimensions. The expectation (O) of local observables O, 
(O) = Tr (O exp{—/JiJD/Z, is then a ratio of contrac¬ 
tions of two TNs. This efficient representation serves as 
the starting point for a coarse-graining process for com¬ 
puting physical observables from the tensor. The con¬ 
traction of a TN is in general a hard problem, and nu¬ 
merical approximations are in order"*^. In recent years a 


number of schemes have been proposed to approximate 
Z in an efficient manner; see e.g. Refs.®^’*® and refer¬ 
ences therein. The RG process consists of a decimation 
of the lattice at each step and such decimation process 
is iterated, which gives rise to efficient and accurate pre¬ 
dictions of thermodynamic quantities, even close to the 
phase transitions. 

Let us take for example the Ising model with nearest- 
neighbor interaction and a local field, 

H = - —{si + Sj)] (10) 

d,.-) 

where {i,j) represents the nearest-neighbor pair i and j, 
Ub is the number of neighbors, nt = 2d for square (d = 2) 
and cubic (d = 3) lattices. Its partition function Z can 
be expressed as a tensor network; the explicit description 
of the tensors forming the network are obtained straight¬ 
forwardly from the Hamiltonian, 

Z = Tr exp(—/3id) = Tr exp{—(11) 

can be decomposed as the tensor trace of the product of 
identical tensors T, 

Z = tTT:\YTT ...T, (12) 

where tTr denotes the tensor trace, i.e., summing over all 
degrees of freedom, and each tensor T (lying on each of 
the edges of, e.g., the square lattice, see Fig. 1) has the 
expression 

/exp{/3-bd/d} exp{-/3} \ , . 

exp{—/3} exp{/3 — d/d} J * ' 

The product operation transforms the tensor network 
into a number by contracting i.e. summing over all the 
corresponding degrees of freedom of neighboring sites. 
Each tensor T has rank 2 and is located on an edge of the 
square lattice. For convenience of an RG iteration that 
directly preserves the square and cubic lattice structure 
of the tensors, we shall instead construct new tensors A 
that are located on each site and the trace of them is to 
sum over degrees of freedom on edges. To do this, we 
first use a simple decomposition (such as singular-value 
decomposition) of each tensor T as 

(14) 

and then we form on each vertex a new tensor A from 
the contraction of four tensors (chosen from U or V) in 
the square lattice 

Au,d,l,r = U^,uV,,dV,,lU,,r- (15) 

i 

(For the simple cubic lattice, A will be a rank-6 ten¬ 
sor obtained from six U or V tensors.) Therefore, the 
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FIG. 1. (Color online) The free energy of a classical lattice 
system is expressed as a Tensor Network, and evaluated using 
a RG iterative process, (a) The Tensor T of rank 2 is obtained 
from the Hamiltonian (see Eq. 13 in the text). By means 
of a decomposition (b), tensors U and V are obtained, and 
combined (c) to form tensor A. An RG step over the square 
lattice consists in the coarse grain of a set of tensors into a 
single site tensor, /.e. combining two tensors horizontally 
(d), a new tensor A is created (e). After truncation, the new 
tensor (f) is used again in an iterative process. 



FIG. 2. (Color online) The HOTRG algorithm proceeds at 
each step combining 2 tensors along a preferred direction, and 
extracting unitaries on the combined indices. These unitaries 
W are truncated and inserted in order to reduce the dimen¬ 
sion of the auxiliary indices. This prevents an explosion of 
the number of parameters required, while keeping accurate 
approximations to Z. Similar operations are required either 
for 2D lattices (top, with tensors of rank 4) and for 3D lattices 
(bottom, tensors of rank 6). 


partition function of the Ising model on a square lattice 
involves a single tensor A repeated in the infinite lattice, 
and is an exact representation of the partition function 
of the system, which is then evaluated via an RG itera¬ 
tive process to be described below. The basic structure 
of an RG process to compute Z is a partition of the lat¬ 
tice into small plaquettes. At each step, plaquettes are 
coarse-grained to obtain an effective lattice of smaller 
size. In the limit of many iterations of this process, pla¬ 
quettes become invariant upon further renormalization 
step. At this effective thermodynamic limit we can eas¬ 
ily evaluate the free energy per site or local observables 
using a single plaquette. While the general picture of 
RG methods is essentially equivalent but differs in the 
specihc algorithmic implementations, in this paper we 
employ the higher-order tensor RG (HOTRG)^®, which 
among other schemes has shown accurate evaluation of 
physical observables. The basic ingredient of the algo¬ 
rithm is presented in Fig. 2, where a pair of tensors A 
is coarse-grained into a single tensor A. In practice, the 
coarse-graining procedure will be done iteratively for the 
horizontal direction followed by the vertical direction in 
the 2D square lattice, and similarly for the three direc¬ 
tions in the 3D simple cubic lattice. 

From the RG flow one directly obtains approximations 
to the free energy per site''’° 


/ 




(16) 


Before starting the RG, we have a TN representing the 
Ising model in a square lattice with initial equal tensors 
A(0); as we perform the RG contraction numerically, we 
keep the tensors under machine representation factoriz¬ 
ing at each step as follows, 

( 17 ) 


From this the RG process effectively shrinks the lattice 
sites by half and produces a sequence of new tensors 
A(i), and each of these tensors is also factorized 

accordingly to bound the numerical values similar to the 
above Eq. (17). Let us define 

G(°)=ln|ao| ( 18 ) 

and the corresponding The partition function in 

the TN picture then reads 

Z=[A(o)]^ = e'^G™[A(°)]^, (19) 

where we have defined the notation [A]-^ to indicate a 
contraction of tensor network described by A over N 
sites. By applying an RG step over A^^^ we have 

Z = (20) 

where A^^\ after a single RG step, spans a smaller lattice 

or more precisely half the lattice with N/2 sites. This 
process can be iterated and the free energy per site can be 
written as a function of all the prefactors as follows, 

( 21 ) 

k=0 

where n the total number of RG steps that has been 
carried, and the second term vanishes exponentially for 
large n. Notice that without truncation the dimension 
of each index of the tensor A will increase exponentially 
with n, and the common practice is to limit the maximal 
dimension to some number denoted by ■ One can 

increase D^ut to see whether the computed observables 
converge. By the above procedure, we thus obtain the 
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free energy per site (at any complex field h value) solely 
from the evaluation of the prefactors along the RG 
flow. In our calculations, the estimation of the free energy 
converges after only a few RG steps (typically n ~ 20). 
We note that previous application of the tensor network 
methods for / has been focused on real h values and high 
precision for observables has been obtained'^^’^"*. 

Would the TN methods work for complex field val¬ 
ues? We use the above method to perform a calculation 
of the free energy per site / on the complex plane de¬ 
fined by z = exp{—2/?/i}. The result is shown in Fig. 3 
for the Ising model in a 2D square lattice, for illustra¬ 
tion, at two different temperatures. For any tempera¬ 
ture, one always observes the minimum of / at exactly 
the unit circle, a verification of the Lee-Yang theorem. 
For T < Tc, the minimum of / is located exactly along 
the unit circle, with a uniform density for any value of 
9 with z = r exp(i0) at r = 1. By increasing T above 
Tc = 2/ log(l -|- v^), the values of / around the positive 
real axis start to increase, and no discontinuity is ob¬ 
served near 0 = 0 at the unit circle; observe that the dark 
region recedes from the real axis close to unity. From the 
point of view of the magnetization (to be discussed in 
the next section), this translates into the disappearance 
of phase transition at this high temperature. These re¬ 
sults demonstrate how tensor network methods can be 
useful in the study of partition function zeros, via the 
free energy / on the complex plane, obtained after an 
efficient RG contraction process. In the next section we 
shall show that the computation of magnetization leads 
to a more precise probe of the location of zeros and their 
density. 


IV. MAGNETIZATION AND DENSITY OF 
YANG-LEE ZEROS IN THE ISING MODEL 

While one can obtain physical quantities directly from 
the derivatives of the free energy, the RG algorithms also 
provide a direct approach to compute expectation values 
of local observables such as the energy and magnetiza¬ 
tion. In order to compute the magnetization we define a 
new tensor B 

i 

where mi = -l-I and m_i = —1, accounting for the local 
magnetization on a single site. Thus, calculations of the 
single site magnetization imply evaluation of the expres¬ 
sion 

m = ^Tr(sie“^^) = (23) 

involving two contractions of a TN: one contraction com¬ 
putes the norm Z (using exactly the tensors defined for 
the calculation of the free energy), while the second con¬ 
traction differs from the first one by only the single site 
tensor defined in Eq. (22) at only one site. We note 
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FIG. 3. (Color online) Free energy per site / in the com¬ 
plex plane defined hy z = exp (— 2/i/3) for 13 = Pc. (top) and 
P = Pc/‘i- (bottom), as obtained from Tensor Network calcu¬ 
lations with Dcvit = 10. Darker regions around the unit circle 
correspond to minimum values of /. For j3 < Pc this mini¬ 
mum does not contact the real axis at z = 1. The central 
region possesses very large values of free-energy density and 
is removed so to enhance contrast around the unit circle. 


that due to the complex held, the magnetization m is 
not necessary real and nor restricted in the range [—1,1]. 
Along the RG process, (M) is balanced by a prefactor 
jafej, obtained from the RG how from Z, and that keeps 
the numerical process bounded at each step. We use the 
ratio Eq. (23) along the RG procedure to ensure conver¬ 
gence, which normally takes about 50 steps. Performing 
two contractions after the RG has converged, we obtain 
a direct measurement of the magnetization in the state 
represented by the TN. This allows a complete study of 
the magnetization for all complex values of the magnetic 
held h. 

Eor the Ising model on the 2D square lattice, only two 
exact solutions of m are known in the complex plane'^^: 
(i) the seminal solution at z = 1, 


m(z = 1) 


1 + x'^ 

(l-a;2)2 


(1 — + x^) 2 ) 


1 

4 


(24) 


and (ii) on the opposite side of the unit circle, at z = —1, 

"(1 -I- x^)^ 


i(z = -1) = 


-(l + 6x^ 


1 — x^ 


(25) 
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FIG. 4. (Color online) Magnetization as a function of a: = 
exp (—2/3) at the points z = 1 and z = —1 for the 2D Ising 
model. Solid lines are obtained from the analytical expres¬ 
sions in Eqs. (24) and (25), and marks are obtained using 
Tensor Networks with Dcut = 15. While the precision de¬ 
creases around the transition point, it remains below 10“^ for 
this choice of Dcut (see inset). 

with X = exp (—2/3). We use these exact analytic results 
to benchmark the performance of the RG contraction. In 
Fig. 4 we present calculations of m vs. x at points z = \ 
and —1, and from comparison with the exact results, we 
find that the error of m (as evaluated for different Dcut) 
is bounded below 10“^, even close to the transition point. 
This clearly demonstrates that TN methods can provide 
an accurate picture of the magnetization for complex val¬ 
ues of the magnetic field. We remark that close to a 
phase transition, numerical accuracy can be improved 
by increasing the bond dimension {Dent)- Furthermore, 
more sophisticated RG algorithms can be employed^®’®^. 
Eqs. (24) and (25) already provide useful information re¬ 
garding the properties of the model and the partition 
function zeros. According to Eq. (24), at z = 1 (or equiv¬ 
alently h = 0) there is a critical value of x (or equivalently 
ksTc = ^1 Pc) beyond which no magnetization is present, 
showing a phase transition. This transition for the 2D 
Ising model appears at Pc = In (I -|- •\/2)/2. At 0 = tt, 
however, m always increases with the temperature (with 
a value becoming larger than unity increasingly), and a 
divergence builds up here in the local magnetization at 
very large temperature T. This buildup of divergence 
also implies the buildup of partition function zeros; see 
Eq. (27) and further discussions below. 

In Fig. 5 we plot for the 2D Ising model the magneti¬ 
zation in the complex z plane. It shows a discontinuity 
only at exactly the unit circle. For the magnetization 
at two different temperatures we can observe how the 
temperature affects this discontinuity. For T < Tc this 
discontinuity appears all along the unit circle, even close 
to the positive real axis. However, at a higher temper¬ 
ature T > Tci the discontinuity around 0 = 0 vanishes. 



FIG. 5. Magnetization per site in the complex plane defined 
by z = exp{—2/3/i}, as obtained from the TRG with Dcut = 8, 
for the Ising model on the 2D square lattice. The temperature 
is fixed at P = Pc (top) and P = Pcjl (bottom). A disconti¬ 
nuity in m appears only at points on top the unit circle, and 
for a range of angles 0 that changes with the temperature. 


and a smooth region emerges around there close to the 
positive real axis. 

The RG method presented above can be extended to 
the study of higher dimensional lattices. Using a tensor 
network structure resembling that of the spin lattice, the 
connectivity of each component increases and the rank 
of the tensors becomes larger at higher lattice dimen¬ 
sions. Thus, the computational cost associated to the 
RG process (which depends directly on the rank and di¬ 
mension of the tensors) is larger, and we can only reach 
smaller bond dimensions (and hence less precision) using 
the same resources, compared to 2D. For the Ising model 
in a simple 3D cubic lattice the tensors have rank 6, and 
are obtained in a similar way to the 2D case, resulting in 
the tensor 
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Using the HOTRG procedure combining two tensors of 
rank 6 at each step, the 3D lattice gets contracted to a 
single tensor. The magnetization in the complex plane is 
obtained also similarly as in 2D and is plotted in Fig. 6. 



FIG. 6. Following Fig. 5, magnetization per site m in the 
complex z plane for the Ising model on the 3D simple cubic 
lattice. Two temperatures are fixed, j3 = pc (top, Dcut = 
5) and /3 = /3c/2 (bottom, Dcut = 9). Above the critical 
temperature, the discontinuity of m turns into a smooth slope 
around the positive real axis. 

How do we then obtain the density of the zeros? The 
Lee-Yang theorem states that for the ferromagnetic Ising 
models, the partition function zeros lie on a unit circle in 
the complex z plane. The density of zeros g{0), for any 
9 along the unit circle re*® with r = 1 is related to the 
discontinuity of the magnetization'^’^^, 

lim Re{m) — lim Re{m) = —47rg(0). (27) 

r—>-1+ r—¥l~~ 

From our calculations we observe how the magnetiza¬ 
tion shows a discontinuity at the unit circle, and how 
this discontinuity changes with the temperature. For 
the Ising model, we observe that lim^_,.]^+ Re{m) = 



e 


FIG. 7. (Golor online) Magnetization m{r —>■ 1~) for the 
2D Ising model using the TN TRG with bond dimension 
Dcut = 45. Via the relation between the magnetzation and 
the density (28), the main features of the density of Yang-Lee 
zeros are obtained in a range of temperatures. For T <Tc, a. 
uniform distribution of zeros is obtained along the unit circle. 
At T = Tc, a drop in the density near the positive real axis 
marks the presence of a phase transition. For T > Tc, a gap 
around the real axis appears and no zeros are found up to a 
critical value Be- 

— limr_,,]^- Re(m) and hence 

lim Re{m) = 2TTg{9). (28) 

r—>-1“ 

Using Eq. (27) or Eq. (28) we can thus directly relate 
our calculation of m to the density of zeros along the 
unit circle. Confirming the Lee-Yang theorem, our re¬ 
sults show that the density g{9) is zero around the real 
axis for T > T^, and thus forms a gap of zero density 
for 6 < Be, where 9^ the smallest value of 9 below which 
no zeros occur. As the temperature increases, the zeros 
move towards 0 = tt, and the density there keeps in¬ 
creasing with temperature, which is what was revealed 
by Eq. (25). At precisely T = Tc the density at the point 
z = 1 drops to zero, indicating the temperature threshold 
of the existence of phase transitions. We plot in Fig. 7 
the density g{9) for three different temperatures around 
Tc- The density of zeros is flat below Tc, and exactly 
at T = Tc we observe a drop in the density around the 
real axis. An entirely different behavior for the zeros is 
observed at T > Tc. 

At the critical temperature T = Tc the relation be¬ 
tween the magnetization and the external field implies 
the following scaling relation of the density of the zeros, 

g{9) ~ |0|® (29) 

where S is the magnetization-field critical exponent. In 
Fig. 8 we plot the density at T = Tc, where we observe 
the detailed drop of the density at exactly z = 1, and the 
inset shows the estimation of the exponent. From these 
calculations we obtain an estimation of 5 = 15.0(2) in 






























































FIG. 8. (Color online) Magnetization along the unit circle 
(r —>■ 1~) for different values of 9 at exactly /3 = /3c. Via the 
relation between the magnetization and the density (28), the 
inset shows the relation log(m) vs log(^) following the relation 
g ^ 1 ^ 1 ^, with a value of 6{Dcut = 30) = 15.0(2) (for the Ising 
model in 2D, S = 15). 

agreement with the 2D Ising magnetic exponent 6 = 15^. 
Proceeding in a similar way for the higher dimensional 
lattice, we use in our calculations an estimated critical 
temperature = 4.51154®^“®® for the 3D Ising model 
in a simple cubic lattice. Plotting the density at this 
temperature (see Fig. 9) we obtain the critical exponent 
S = 4.8(3) from the relation g ^ |0|'^ (for a 3D Ising 
model in a simple cubic lattice i5 = 4.789(2)). 


V. EDGE SINGULARITIES AT T > Tc 

The structure of the density of zeros changes dramat¬ 
ically at temperatures above as a, gap with no zeros 
on the unit circle opens up around 6 = 0, and it increases 
with the temperature. The density of the Yang-Lee ze¬ 
ros vanishes for 6 < 6e and becomes non-zero starting at 
an edge value The value 6^ can be directly deter¬ 

mined from the density obtained using the TN method, 
and has strong dependence with the temperature. From 
the magnetization we obtain the value of df,{T), depicted 
in Fig. 10 for 2D square and 3D simple cubic lattices, 
as a function of the temperature T normalized to their 
respective Tc- While not exactly equal, these two curves 
show similar features as they both increase rapidly after 
Tc and then have a slower progression at larger temper¬ 
atures. 

In addition to 9e, there is also a global movement of 
zeros away from 9 = 0 towards the opposite side of the 
unit circle 0 = tt as the temperature increases. In Fig. 11 
we plot the density obtained around the unit circle as 
a function of 9 for a 2D square lattice. Eventually at a 
sufficiently large temperature, the density will accumu¬ 
late mostly at 9 = TT and diverge as T —>■ oo. This is 



e 


FIG. 9. (Golor online) For an ising model in a cubic lattice, 
the magnetization and hence the density of Yang-Lee zeros 
at the estimated critical temperature Tc = 4.51154 along the 
unit circle in the complex plane is plotted using a bond dimen¬ 
sion Dcut = 14. Via the relation between the magnetization 
and the density (28), it is observed that this density is zero at 
the real axis 0 = 0, and increases along the unit circle. The 
inset shows the calculation of the critical exponent S = 4.8(3) 
from the relation g ~ |0|'* (for a 3D Ising model in a cubic 
lattice S = 4.789(2)). 


also verified by the diverging behavior of m at z = — 1 in 
Eq. (25). 

The most striking feature in the behavior of the zeros 
is the edge singularity in 2D: as 9 approaches 9^, from 
above, the density of zeros becomes diverging®®. For in¬ 
creasing values of T, the density on the unit circle evolves 
as illustrated in Fig. 11. This singularity at 9e has been 
identified as a critical pointand as a non-unitary re¬ 
alization of conformal symmetry^®. It is characterized as 
follows, 

g{9)^{9-9cY, for0>6»e, (30) 

where 9e is the position of the divergence at each temper¬ 
ature. One expects a reduction of accuracy near a critical 
point, especially for this edge-singularity critical point, 
and we only obtain an estimated value of a = —0.1(1), 
consistent with the value from conformal field consider¬ 
ations (7 = —1/6*^®. A different picture appears in the 
study for the 3D simple cubic lattice. There is no di¬ 
vergence of the density and a is positive®'; see Fig. 6b. 
However, due to small D^ut we can use and the noisy 
data from the density, we do not obtain a good estimate 
of <7. 


VI. YANG-LEE ZEROS IN THE POTTS MODEL 
WITH A COMPLEX EIELD 

In the preceding sections we have presented a method 
to obtain properties of the partition function zeros for 
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FIG. 10. (Color online) Location of the partition function 
zero closest to the positive real axis, as measured by the an¬ 
gle 6s. Results for the 2D square lattice (using Dcut = 20) 
and for a 3D simple cubic lattice (using Dcut ~ 10) are shown, 
normalized by their respective Tc- The location of the edge 
singularity 9c is shifted towards 0 —>■ tt for increasing temper¬ 
atures. 
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FIG. 11. (Color online) Density of zeros for the 2D Ising 
model at different temperatures above the critical value Tc, 
as computed with TN and Dcut = 45. A gap of zeros appears 
around the real axis for a range of 9, e.g., up to a value of 
0.8 for /3 = Pc/3. This gap extends to 9c, the edge singularity 
with a high density of zeros. The value of 6c moves towards 
9 = TT as the temperature is increased. Inset: At /3 = 0.2, 
estimation of the exponent fi = —0.1(1). 


lattice models, focusing on the Ising model in two and 
three dimensions. This approach can indeed be applied 
to any model for which we have a TN description, and 
especially in classical models where the tensors are read¬ 
ily obtained directly from the Hamiltonian expression, as 
explained in the Ising models above. For illustration, we 



FIG. 12. (Color online) The location of the zeros for the 
2D g-state Potts model using polar coordinates, as obtained 
with Dcut ~ 20. At the corresponding critical temperature 
/3c = log(l -f ^/q) for each q, the locus of zeros lies outside the 
unit circle. The location of zeros is shown for the Potts model 
with g = 3 ... 6, where the distance changes with 9 showing a 
maximum value at 9 — n for all g. 

shall consider the g-state Potts models on the square and 
simple cubic lattices. 

The tensor network description of the g-state Potts 
model is obtained from 

H = '^[1 - 6{ai,<Tj)] - h'^S{a^,t) (31) 

ihj) * 

where cr = 0,... g — 1 and t G [0, g — 1] but can be chosen 
to be 0 for simplicity. The transition temperature of the 
g-Potts model was found to be at (3c = log(l -I- 
Following similar considerations as in previous sections, 
we obtain a tensor network of an initial bond dimension g 
on the square lattice (as well as the simple cubic lattice), 
where tensors are located at the vertices of the lattice. 
The TRG method is directly applied to this tensor struc¬ 
ture. We compute the magnetization 

m = ^Tr(Me-'^^), (32) 

Zj 

where Z = Tr exp(—/3H) and 

M = ^^,5(cri,t). (33) 

i 

Similar to the explanation for the magnetization in the 
Ising model (23), the compuation of the magnetization 
in the Potts model is also a ratio between two tensor- 
network contractions. The density of zeros is also ob¬ 
tained by searching for discontinuity in the real part of 
the magnetization on the complex plane. 

The exact location of Yang-Lee zeros of the Potts 
model has attracted attention as we lack the equivalence 
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FIG. 13. (Color online) Distance to the unit circle Ar in the 
complex plane at the point 6 = n for the 2D g-state Potts 
model. A bond dimension Dcut = 20 is used, and distances 
are shown for a number of temperatures /3. Increasing the 
value of q increases the distance Ar, but decreases at higher 
temperatures. However, this plot suggests that only in the 
limit T = 0 the zeros reach the unit circle. 


of the unit circle theorem for this model. Their loca¬ 
tion has been mostly estimated from finite-size extrapo¬ 
lation^®’®®. An interesting feature of these findings is that 
for the Potts model for g > 2 the zeros were believed to 
lie outside the unit circle. Upon increasing the tempera¬ 
ture, the zeros move further away from the unit circle as 
the gap around real positive axis opens. However, these 
previous results based on finite-size extrapolation were 
questioned"^®, and calculations directly in the thermody¬ 
namic limit were not available. 

Our approach explained earlier enables us to probe di¬ 
rectly the thermodynamic limit and to resolve the debate. 
In Fig. 12 we plot the location of zeros (obtained via the 
discontinuity in the magnetization) for the Potts model 
at /3 = /3c, for different values of q. We clearly observe 
how the locus of zeros is outside the unit circle, at a vari¬ 
able distance dependent on 0 and q. We have verified that 
our results have little dependence on the bond dimension 
(for Dcut ^ 20), and this shows that the RG schemes, 
such as the HOTRG, do provide efficient method to ac¬ 
cess the thermodynamic limit. These results are in good 
agreement with those in^®, providing a complete picture 
of the locus of zeros at any 0 in a precise way. As a con¬ 
sequence of the large though finite correlation length at 
g > 4 (of a few thousand sites®®), we use a larger bond 
dimension to achieve similar precision level for any value 
of g. As we observe in Fig. 12, the farthest zero is located 
at 0 = TT in agreement with previous finite-size study^®. 
In Fig. 13 we show the movement of the zeros located at 
this point as measured by Ar = |r — 1|, as we lower the 
temperature (i.e., increase /3). At the limit T = 0 the 
zeros lie on the unit circle, but this limit is reached only 



FIG. 14. (Color online) The location of the zeros for the 
g = 3, 4 Potts model using polar axis in the 3D simple cubic 
lattice, as obtained up to Dcut = 12. At their respective 
critical temperatures®® the locus of zeros lie clearly outside 
the unit circle. 

exponentially slowly with the inverse temperature, at a 
similar rate (Ar ~ for any g and other values 

of 0 (not shown here). 

In three dimensions, previous finite-size study was lim¬ 
ited to very small system sizes such as 3 x 3 x 3®®, but 
the zeros do not lie on the unit circle. What is their fate 
in the thermodynamic limit? Again we use TN methods 
to directly probe the zeros in this limit. In Fig. 14 we 
show the locus of zeros for g = 3,4 Potts model in the 
3D simple cubic lattice. These zeros are computed at the 
critical temperature T = Tc from Ref.®®, and we clearly 
see that they are located outside the unit circle in the 
complex z plane (except at 0 = 0). 


VII. CONCLUSIONS 

We have employed the tensor renormalization meth¬ 
ods to investigate the properties of the free energy and 
the distribution of the partition function zeros (i.e. the 
Yang-Lee zeros) in the complex field (analagous to the 
fugacity) plane. From the position and the density of 
Yang-Lee zeros one can determine important properties 
of the phase transitions of spin systems. While the loca¬ 
tion of the zeros was rigorously established by Lee and 
Yang for the Ising model, their distribution has not been 
established precisely. We have demonstrated that the 
tensor network methods provide useful tools to access 
such information on the complex plane. 

We have presented results for the density of zeros along 
the unit circle in the plane of complex z = exp{—2/3/i} 
in both 2D and 3D Ising models, showing different char¬ 
acteristics of the density at different temperatures com¬ 
pared to the critical temperature. In particular from 
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the distribution of zeros at we have extracted the 
magnetization-field critical exponent. At higher temper¬ 
atures we have determined how the singularity edge 9^ 
moves with the temperature and estimated the singular¬ 
ity exponent. All the results are in good agreement with 
those from other techniques. 

Going beyond the Ising models, we have also exam¬ 
ined the g-state (with q > 2) Potts model in both two 
and three dimensions, where fewer analytic results were 
known. We found that in the thermodynamic limit the 
Yang-Lee zeros from these small system sizes do not lie 
on the unit circle except at the zero temperature, and 


the approach to the unit circle from high temperatures 
is exponentially close as the inverse temperature. This 
resolves a previous debate about whether the conclusion 
that the zeros are not on the unit circle is indeed correct 
in the thermodynamic limit or simply due to the finite- 
size effect. Possible future directions include the appli¬ 
cation and generalization of the approach here to other 
models for probing both Yang-Lee and Fisher zeros. 
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